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ABSTRACT 



We have calculated several representative models of vertical structure of an 
accretion disk around a supcrmassivc Kerr black hole. The interaction of radiation 
and matter is treated self-consistently, taking into account departures from LTE 
for calculating both the disk structure and the radiation field. The structural 
equations are described in detail, and various approximations are discussed. We have 
demonstrated that departures from LTE are very important for determining the disk 
structure, even at the midplane, as well as the emergent radiation, particularly for 
hot and electron-scattering-dominated disks. We have shown that at least for the 
disk parameters studied in this paper, NLTE effects tend to reduce the value of the 
Lyman jump with respect to the LTE predictions, regardless whether LTE predicts an 
emission or absorption jump. Wc have studied the effects of various values of viscosity 
on the model structure and predicted spectral energy distribution. The viscosity is 
parameterized through a parameter ao which describes the vertically-averaged viscous 
stress, two power-law exponents Co and Ci, and the division point md between these 
two forms. The disk structure and emergent radiation is sensitive mainly to the values 
of ao, while the other parameters influence the disk structure to a much lesser extent. 
However, although the detailed shape of the predicted spectrum is sensitive to adopted 
value of ao, the overall appearance of the spectrum is quite similar. 



Subject headings: accretion, accretion disks — galaxies: active — galaxies: nuclei — 
radiative transfer 
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1. INTRODUCTION 

Accretion disks around massive black holes have long been the most popular candidates 
for providing the ultraviolet and soft X-ray flux observed in Active Galactic Nuclei (AGN). 
Observational evidence is mainly based on the 'big blue bumps' seen in the UV (e.g. Shields 
1978; Malkan & Sargent 1982). However, despite its attractiveness this model faces a number of 
problems when confronted with multi-wavelength observations. The current situation has been 
recently summarized by Koratkar (1998) from the observational point of view, and by Blaes (1998) 
from the theoretical viewpoint. The most pressing problems of the theoretical models are a near 
absence of observed Lyman discontinuity (first pointed out by Antonucci, Kinney, & Ford 1989), 
the UV/EUV continuum polarization (e.g., Antonucci 1992, Blaes 1998 and references therein), 
the overall continuum spectral energy distribution (e.g. Laor 1990), and phased optical/UV 
variability (e.g., Alloin et al 1985). 

What is clearly needed is a self-consistent model which would explain all the observed 
features. Such a goal is still far away, but we feel that the first step towards it is to answer the 
fundamental question: In view of all current problems, is the accretion disk paradigm still a viable 
model? In other words, are the observed phenomena truly inconsistent with the accretion disk 
picture in general, or is the lack of agreement rather a result of inaccuracies or even inconsistencies 
in the computed models or in the current modeling techniques? 

Therefore, we have embarked on a systematic study of these questions. We recognize that 
there are many possible sources of inconsistencies in the AGN accretion disk modeling. There are 
essentially two types of problems. First, there arc basic physical uncertainties. Among them, the 
most important is our ignorance of the basic physics of viscous energy dissipation in the AGN 
disks, so that we are left with a necessity to employ certain ad hoc parameters. Although this is 
certainly a viable approach when nothing else is currently available, we should bear in mind that 
corresponding models will lack predictive power — we may explain what we see, but we cannot 
predict it. In any case, when adopting a model based on some chosen set of parameters, one 
at least has to study carefully the influence of these parameters on the computed model. Other 
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fundamental physical problems are uncertainties of the effects of other structures forming the 
AGN complex onto the accretion disk; for instance an irradiation of the disk from external sources. 
Again, in the absence of any detailed theory, we are usually left with a necessity to parameterize 
these effects. 

Second type of problems concerns the degree of approximation used in the actual modeling 
procedure. A particularly important class of such approximations deals with the description of 
interaction of radiation and matter in the disk. It should be emphasized that AGN disks, like 
atmospheres of hot stars, are typical examples of a medium where radiation is not only a probe of 
the physical state, but in fact a crucial constituent. In other words, radiation not only carries an 
information about the medium, it in fact determines its structure. Consequently, a treatment of 
this interaction is in a sense the very gist of the problem. We should therefore study the influence 
of various approximations, such as the degree of equilibrium assumed, or, more specifically, the 
extent to which the local thermodynamic equilibrium (LTE) applies, and the completeness of 
considered opacity and emissivity sources on the computed model, etc. 

In the previous paper (Hubeny & Hubeny 1997 - hereafter called Paper I), we have presented 
some representative self-consistent, non-LTE models of the vertical structure of AGN accretion 
disks. The basic aim of that study was to investigate the differences in the predicted spectrum 
between this and previous approaches (Sun & Malkan 1989; Laor Sz Netzer, 1989; Ross, Fabian, 
& Mineshige 1992; Wehrse et al. 1993; Storzcr and Hauschildt 1994; Coleman 1994; Shields k 
Coleman 1994; Blaes & Agol 1996; Dorrer et al. 1996). The emphasis was to clarify the role of 
departures from LTE and to study the effects of simplifications of the hydrostatic equilibrium 
equation based on assuming a depth-independent vertical gravity acceleration. 

In the present paper, we will consider models of vertical structure of AGN disks in detail. 
In particular, we will study the dependence of computed vertical structure and corresponding 
emergent spectrum on the adopted value of viscosity and on the degree of sophistication of the 
modeling procedure. In order to better emphasize the observable consequences of self-consistent, 
non-LTE models, we present the predicted spectra for a few representative points on the disk. 
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Complete spectra which are obtained by integrating the local spectra over the disk surface, taking 
into account general relativistic photon transfer functions (Cunningham 1975; Speith, Riffert, & 
Ruder 1995; Agol 1997; Agol, Hubeny, & Blaes 1998), will be considered in subsequent papers of 
this series. 

Our aim here is not to construct a model to be used for comparison with actual observations. 
Instead, we intend to study the sensitivity of computed models on the degree of approximations 
used in the modeling procedure. Therefore, we will first give a detailed overview of the structural 
equations and the modeling procedures. This is meant to provide a firm framework on which 
our approach is based, which in turn is useful to assessing possible systematic effects within our 
models. 

2. BASIC EQUATIONS AND APPROXIMATIONS 

2.1. Relativistic Disk Structure 

We assume a steady-state, thin, rotationally symmetric disk. We define the thin disk by the 
two conditions, namely a) the disk height, h{R), at radial distance R is much smaller than R, i.e., 
h{R) <C R, and b) all the components of the energy flux vector are negligible to that in the vertical 
direction, i.e., the direction perpendicular to the disk plane. Vertical distance from the disk plane 
is denoted z. We also assume the standard disk model, in which the azimuthal velocity is much 
larger than the radial velocity, and the vertical velocity component is neglected altogether. 

The equations for relativistic disk structure were derived by Novikov &: Thorne (1973), Page 
& Thorne (1974), Eardley & Lightman (1975) - who have corrected a previously incorrect term in 
the vertical pressure balance, and recently by Riffert &; Herold (1996) whose results we use here. 

The four basic equations describing the disk structure are: 

i) Vertical hydrostatic equilibrium, 

dP GM C 

5J = -^"^5' 
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where P is the total pressure, p is the mass density, G is the gravitational constant, M is the 
black hole mass, and B and C (together with A and D used later) are the so-called relativistic 
corrections - see below. We note that Abramowicz, Lanza, &: Pcrcival (1997) have recently 
rederived the vertical hydrostatic equilibrium equation, and showed that previous treatments 
yielded some unphysical singularities close to the last stable orbit. 

ii) Energy balance equation, 



dF,_3 GMA 

where is the 2;-component of the energy flux, and is the sheer stress (also called the viscous 

stress). 

iii) Azimuthal momentum balance (written using the equation of continuity, and with the boundary 
condition t^r = at the innermost stable orbit) , 



1^ M CM D 



h 



where M is the mass accretion rate. 

iv) Finally, the equation describing the source of viscous stress, 



3 CM A 

V=2^y^5, (4) 
which specifies the viscous stress in terms of velocity gradients and the sheer viscosity, rj. 
The relativistic corrections are given by 

1 r x2 - 6a; + 8a - Sa^ , 
D = — ^ / ^ ^ ^ dx , (8) 
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where r is radius of the annulus expressed in units of the gravitational radius, = GMjc? ^ i.e., 
r = R/ {GM / c^)] and the specific angular momentum a/M is expressed in units of G/c; c being 
the speed of light, and rj the radius of the innermost stable orbit. 

Equations (Q) - (^) are general equations which describe the structure of the disk. To make 
them applicable to real disks, one has to specify, in addition, a) the nature of the energy flux 
(i.e. how the energy generated by viscous sheer is transported), and b) the nature (and numerical 
value) of the sheer viscosity. 

The first item is straightforward. It is usually assumed that the energy is transported solely 
by radiation (plus possibly in part by convection in convectively unstable layers). In this case 
is the total radiation flux. From Eq. (|2|) it follows that 

^ , , , 3 IgM a ^ 

because the flux vanishes at the disk plane, Fz{z = {)) = 0, since the disk is symmetric about the 
disk midplane. It is customary to express the total energy flux at the disk surface through the 
effective temperature. As follows from Eqs. (^) and (^), 

/ , N rr.A 3 GMM D , , 

where a is the Stefan-Boltzmann constant. The effective temperature has thus the meaning 
of the temperature for which the corresponding black-body radiation has the total radiation 
energy (integrated over all frequencies) equal to the total energy generated in the column of unit 
cross-section in the disk. This explains why the early approaches used the black-body radiation 
for describing the disk radiation; however, we stress that the radiation emergent from the disk 
does not have to possess the black-body frequency distribution. 

2.2. Viscosity 

Viscosity is the most uncertain physical quantity of the accretion disk modeling. There is no 
theory that would explain the accretion disk viscosity form first principles, although, at least in 
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the case of accretion disks in the cataclysmic variables, the Balbus-Hawley (1991) instability and 
subsequent detailed numerical MHD simulations (e.g., Stone et al. 1996) represent a breakthrough 
in our understanding of disk viscosity, and offer a great promise for the future. 

In this paper, we adopt the traditional approach, and parameterize the viscosity by means of 
some adjustable parameters. First, we express the sheer viscosity through the kinematic viscosity, 
w, as 

rj = pw . (11) 
We also introduce the vertically averaged kinematic viscosity, 



_ /o wpdz 

w — — 



I rh. I mo 

= — / r]dz = — / wdm, (12) 
m-o Jo mo Jo 



Jq pdz 

where we have introduced the column mass, dm = —pdz, i.e., m{z) is the column mass above the 
height z. The total column density at the disk midplane is denoted by mo, and is related to the 
traditional disk surface density, E, by mo = E/2. 

The most commonly used viscosity parameterization is the so-called a-prescription (Shakura 
&; Sunyaev 1973). In this model, viscosity is thought to be caused by turbulence; the kinematic 
viscosity is then postulated, by a simple dimensional analysis, to be equal to 

^i' = a ^turb ^^turb , (13) 

where Zturb and Uturb ^.re the size and velocity of the largest turbulent cells, respectively, and a is 

an ad hoc proportionality constant. There are several variants of the a-prescription, depending of 

what is taken for Zturb and wturb- Typically, kurh is taken to be equal to h, the disk height, and 

■J^turb equal to the sound speed, Cg = \/ (P/p), or some other turbulent velocity. This is a local 

description. In fact, the original Shakura-Sunyaev a was introduced to describe vertically-averaged 

quatities, and we use this description here as well. We define 

I fh _ 

= r / Ur dz = aoP , (14) 
ri Jo 

where P is the vertically-averaged total pressure. Further, we write 

rh _ _ 

tsr dz = haoP = moao (P/p) , (15) 
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where p = niQ/h is the verticaUy averaged density. 

In this paper, we do not consider P and p to be the model-dependent averages of actual 
pressure and density. Instead, the factor (P/j)) is taken as a known function of the basic 
parameters of the disk and the radius of the annulus, having the value corresponding to the case 
of radiation-pressure-dominated disks, (P/p) = {P/p)iaA- This value is given by (see Sect 3.1 for a 
detailed derivation) 

P _ 3GMM f a, y 

p i?3 [sTTmHcJ BC ^ ' 

where is the electron (Thomson) scattering cross-section, and mn is the mass of hydrogen 
atom. The vertically averaged kinematic viscosity then follows from integrating Eq. ^ over z, 
and using Eqs. (|T2|), (jl|) and (|l|), 



w = 2M^ao\ — 5- — — . (17) 

V \8TTmHcJ AC ^ ' 

The advantage of this parameterization of viscosity is that it yields the total column mass as an 
explicit function of radial distance. Substituting Eqs. (|T^ and ( [T^ ) into (Q), we obtain 



Because of this feature, we can easily use the mass column density, m, as an independent depth 
variable of the problem. This has a significant benefit of greater numerical stability of the solution 
of the structural equations, in particular the radiative transfer and the hydrostatic equilibrium 
equations - see below. 

Alternatively, one may parameterize the vertically-averaged kinematic viscosity through the 
Reynolds number. Re, which is an approach suggested already by Lynden-Bell & Pringle (1974), 



in which case 



^ = ^^, (19) 



MRe BD 
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which fohows directly from Eqs. and (^). The relation between follows from 

Eqs. dig) and (0), 

V f^c / VM/ ao ^2 

The (depth-dependent) viscosity w is allowed to vary as a step-wise power law of the mass 
column density, viz. 

w{m) = wq {m/mo)^° , m > rud , (22) 

w^m) = Wl {m/mo)^^ , m < , (23) 

where ma is the division point. In other words, we allow for a different power-law exponent 
for inner and outer layers. This represents a generalization of an approach we used previously, 
based on a single power-law representation introduced by Ki^fz & Hubeny (1986). The reason 
for choosing a two-step power law is that with a single power law we typically obtain a density 
inversion in the deep layers, while the main reason for introducing the power-law representation of 
viscosity in the first place was to prevent the "thermal catastrophe" of the disk in the low optical 
depth regions where the cooling due to strong resonance lines of light metals is important (see, 
e.g., Shaviv & Wehrse 1986; Hubeny 1990a). On the other hand, recent numerical simulations 
(Stone et al. 1996) indicate that that the viscosity actually increases towards the surface, giving 
rise to a high-temperature region (corona) on the top of a disk. Such models were considered for 
instance by Sincell & Krolik (1997). We plan to study such non-standard models in future papers 
of this series; in this paper we will consider models where the viscosity decreases in the outer 
layers. 

We have thus four independent parameters: exponents Co and C,i, the division point, md, 
and the fraction, /, of energy dissipated in deep layers, m > mj. The coefficients wq and wi 
are derived from the condition on the vertically averaged viscosity, J^" w{m)dm/mo = w, and 
Xnd w{m)dm/mo = fw. We obtain 

fwjCo + 1) 

= 1 - (md/mo)Co+i ' (2^) 

Wl = (^-^)^(^^ + ^) . (25) 
(md/mo)^i+i 
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Generally, w{m) does not have to be continuous at the division point rrid- If we require the 
continuity, then / and are no longer two independent parameters; instead, they are related 
through 



Typically, the deep-layer power law exponent is set to (constant viscosity), while the "surface' 
power law exponent is usually set to a value larger than zero. In §3.3, we will compare this 
treatment to a purely local a-parameterization of viscosity. 



Due to the assumption of thin disk, we may reduce a general 2-D problem of computing disk 
structure to a set of 1-D problems. The disk is divided into a set of axially symmetric concentric 
annuli; each annulus behaves as a one-dimensional radiating slab. The vertical structure of a 
single annulus is computed by solving simultaneously the hydrostatic equilibrium equation, the 
energy balance equation, the radiative transfer equation, and, since we do not generally assume 
LTE, the set of statistical equilibrium equations. Below, we specify these equations in detail. Since 
the state parameters now depend only on the vertical distance z, we replace partial derivatives in 
Eqs. (|l|) and (|2|) by ordinary derivatives. Moreover, we take the column density m as the basic 
depth variable, as it is customary in modeling stellar atmospheres and non-relativistic accretion 
disks. Furthermore, we write down the terms corresponding to radiation (i.e., radiation pressure, 
and radiation flux) explicitly, to stress that we treat the radiation essentially exactly, without 
any simplifying assumptions about radiative transfer. Since the radiation terms in the structural 
equations are written using the radiative transfer equation, we start with it. 

a) Radiative Transfer Equation 

The radiative transfer equation is written in the standard way (e.g. Mihalas 1978), viz 




(26) 



2.3. Equations for the Vertical Structure 



or 




(27) 
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where luip) is the specific intensity of radiation at frequency u, /i being the cosine of the angle 
between direction of propagation and the normal to the disk midplane. The monochromatic optical 
depth is defined through dr^, = —Xydz = {xu/p)dm] and the source function by = r]u/Xv- 
Here, Xu is the total absorption coefficient, Xu = i-^u + ctu, (^u being the scattering coefficient. We 
assume that the only scattering process is the electron (Thomson) scattering, = a^; being 
the Thomson cross-section. Finally, r]i, is the emission coefficient, given by r]y = rf^ + a^Jy, 
where "rf^ is the coefficient of thermal emission. We assume the symmetry condition at the disk 
midplane, /jy(^) = Iy{—fi). The upper boundary condition is = , fJ- < 0, where /~ is 

a prescribed incident radiation. In the present paper we assume, for simplicity, the case of no 
external irradiation, = 0. We return to the problem of non-zero external irradiation in a future 
paper of this series. 

The first two moment equations of the transfer equation read 

dHy 1 



dm p 



i^uJu - Vu) , (28) 



'""^ ^Hy, (29) 



dm p 

where the moments of the specific intensity are defined by 

[J„ H„ K,] = ^ J\l, fi, fi^] /,(^) d^ , (30) 

We stress that in Eq. (28) the scattering terms cancel, so that the equation contains the true 
absorption coefficient k,^, instead of the total absorption coefficient Xu- The radiation flux is given 
by 

POD 

Frad = 47r / Hydiy . (31) 







To solve the radiative transfer equation, we employ the Variable Eddington Factors technique 
(Auer & Mihalas 1970), which consists in introducing the Eddington factor = Ky/Jy, and 
writing the transfer equation (p7|) as 



dT?, 



Jy-Sy. (32) 
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This equation involves only the mean intensity of radiation, Jj^, which is a function of only 
frequency and depth, and not the specific intensity which in addition is a function of angle /x. Such 
an approach is extremely advantageous in methods which solve the global system of structural 
equations iteratively. The Eddington factor is determined by a set of frequency-by-frequency 
formal solutions of the transfer equation, and is held fixed at a subsequent iteration step of the 
global scheme. 

b) Hydrostatic Equilibrium 

The hydrostatic equilibrium equation reads, neglecting self-gravity of the disk and assuming that 
the radial distance from the black hole, R, is much larger than the vertical distance from the 
central plane, z, 

where the depth-dependent vertical gravity acceleration is given by 

t ^ GMC 

5(^) = -Rr:B^- (34) 

The total pressure is given as a sum of the gas pressure and the radiation pressure, 

47r 

P = -Pgas + -Prad = NkT +— du , (35) 



c Jo 

where is the total particle number density, T the temperature, k the Boltzmann constant. The 
upper boundary condition is taken from Hubeny (1990a - Eqs. 4.19-4.20 there). 

c) Energy Balance 

Substituting Eqs. (|) and (|ll]) into (||), and using Eqs. (^) and (|^, we obtain 

9 GM f AX"^ 

A~W\b) P^ = ^'^J^ {r]u- K-uJu)du. (36) 

The energy balance equation may be cast to different form if we do not express the radiation flux 
through the moment equation of the transfer equation, namely 

- wim) . (37) 



dm 4 E? \B 
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Since the radiation flux at the midplane, m = mo, is zero, equation (^7|) may be integrated to yield 



mo 



^ , , 9GM ,,,,,, 

i^radM = --^ f -j / «;(m')dm'. (38) 



m 



Using Eqs. ( p^ ) - pO|), and the definition of effective temperature, Eq. (10), we obtain finally 

i^rad(m) = aTjj [1 - eim)] , (39) 
where the auxiliary function 9 is defined by 

1 Z"™ 

e{m) = - / w(m')dm'. (40) 

wmo Jo 

For any depth dependence of viscosity, is a monotonically increasing function of m, between 
6(0) = 0, and 0{mo) = 1. For a depth-independent viscosity, 9{m) = m/mQ, and for the adopted 
step- wise power law defined by Eqs. (p^) - (p6|), it is given by 



9{m) = {1 - f) (m/md)'^^'^^ , m < rrid , (41) 
af \ (^ f\ I f i'm/mo)<°+^ - (md/mo)^P+^ 

9{m) = [1 - f) + f — — , m>md. 42) 

1 — [md/moj'^o+i 

d) The z-m relation 

Since the hydrostatic equilibrium equation (^) contains the vertical distance z explicitly, we have 
to supply the relation between z and m, which reads simply 

(43) 

am p 

e) Absorption and emission coefficient 
The absorption coefficient is given by 

Xu = J2T.i''i-(9i/9j)nj]a^j{l^)+J2{ni-n*e-''''/''^)a^,{l^) 

i j>i i 

+ ^ n,n^a^^{u, T) (l - 6"'^'^/'=^) + n,a, , (44) 



where Ui is the number density (population) of an energy level i (we number all levels consecutively, 
without notational distinction of the parent atom and ion), n| the corresponding LTE population, 
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gi the statistical weight, and cr(i^) the appropriate cross-section. The four terms represent, 
respectively, the contributions of bound-bound transitions (i.e. spectral lines), bound- free 
transitions (continua), free- free absorption (inverse brehmstrahlung) , and of electron scattering. 
Subscript k denotes the "continuum" , and the ion number density. The negative contributions 
in the first three terms represent the stimulated emission. There is no stimulated emission 
correction for the scattering term, since this contribution cancels with ordinary absorption for 
coherent scattering (for an illuminating discussion, see Shu 1991). 

Analogously, the thermal emission coefficient is given by 

i j>i i 

+ Enen«c7«,(i.,r)e-'^'^/*^^]. (45) 

K 

The three terms again describe the bound-bound, bound-free, and free-free emission processes, 
respectively. 

These equations should be complemented by expressions for the relevant cross-sections, 
definition of LTE populations, and other necessary expressions. 

We did not yet implement the Compton scattering in our codes, but the work on this problem 

is under way, and will be reported in a future paper. Nevertheless, Compton scattering is not 
important for the models considered in this paper, as we shall verify in Sect. 3.1. 

f) Statistical Equilibrium Equation 

It is well known that the LTE approximation breaks down in low-density, radiation-dominated 
media (see, e.g. Mihalas 1978), which are precisely the conditions prevailing in the AGN disks. 
Therefore, we have to adopt a more general treatment, traditionally called non-LTE (or NLTE), 
where the populations of some selected energy levels of some selected atoms/ions are allowed to 
depart from the Boltzmann-Saha distribution. These populations are determined through the 
equations of statistical equilibrium (e.g. Mihalas 1978), viz. 

Hi E {Rij + Cij) = E {Rji + Cji) , (46) 
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where Rij and Cij is the radiative and coUisional rates, respectively, for the transition from level 
i to level j. The l.h.s. of (^) represents the total number of transitions out of level i, while the 
r.h.s. represents the total number of transitions into level i from all other levels. An essential 
numerical complication inherent to this approach is that the radiative rates are given as integrals 
over the radiation intensity, i.e., schematically 

Rij= {hv/A^)aij{v)J^du . (47) 

JO 



2.4. Numerical Method 



The overall system of Eqs. (|32|), (33), (p6|), (^3|), and (|^), together with auxiliary expressions 



(34), ( p!9| ) - (p6[), and the definition expressions for the absorption and emission coefficients, ( ^4[ ) 
and (|45|), form a highly coupled, non-linear set of integro-differential equations. Fortunately, these 
equations are very similar to the equations describing classical NLTE stellar atmospheres (see, 
e.g., Hubeny 1990a, b), where the modeling techniques are highly advanced. We may therefore 
employ to great advantage numerical methods and computer programs designed originally for 
stellar atmospheres. 

We use here the computer program TLUSDISK, which is a derivative of the stellar 
atmosphere program TLUSTY (Hubeny 1988). The program is based on the hybrid complete- 
linearization/accelerated lambda iteration (CL/ALI) method (Hubeny & Lanz 1995). The method 
resembles the traditional complete linearization, however the radiation intensity in most (but not 
necessarily all) frequencies is not linearized; instead it is treated via the ALI scheme (for a review 
of the ALI method, see e.g., Hubeny 1992). A NLTE model of a vertical structure of one annulus 
of a disk is generally computed in several steps. First, an LTE-gray model is constructed, as 
described in Hubeny (1990a). This serves as the starting solution for the subsequent step, an LTE 
model, computed by TLUSDISK. This in turn is used as the starting solution for the next step, a 
NLTE model. 

In the NLTE step, we first calculate the so-called NLTE/C model (i.e., NLTE with continua 
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only) , assuming that all bound-bound transitions are in detailed radiative balance. Finally, in the 
last step, we consider all lines - we denote this model as NLTE/L. The lines influence the disk 
structure both directly - through their contribution to the total heating/cooling rate and to the 
total radiation pressure, as well as indirectly - through their influence on atomic level population 
via equations of statistical equilibrium. Therefore, they influence the heating/cooling rates and 
the radiation pressure in the continuum processes as well. However, one should be aware that 
considering lines is not, strictly speaking, consistent with the 1-D approach adopted here. The 
assumption of horizontally homogeneous rings implies that the Keplerian velocity is assumed 
constant within the ring. This is a good approximation for continua, but may be inaccurate for 
the radiation transfer in spectral lines because a difference of projected Keplerian velocity as small 
as the thermal velocity already shifts the line photon out of the Doppler core. The photon escape 
probability from a disk ring may therefore be quite different compared to the escape probability 
from a plane-parallel, static atmosphere. For a proper treatment of this problems we would have 
to abandon a 1-D modeling and construct a fully 2-D model. However, at present, we intend to 
explore a magnitude of various phenomena and their effect on the disk structure; we feel that for 
such an exploratory model study the 1-D treatment as described above is satisfactory. 

3. SENSITIVITY ANALYSIS 

As a representative case, we take a disk around a Kerr supermassive black hole with 
M = 2 X 10^ Mq, with the maximum stable rotation (the specific angular momentum 
a/M = 0.998). The mass flux is taken M = 1 Mq/jgw. We have calculated a number of vertical 
structure models at various radii; we present here three representative models for r = 2, 11, and 
20. The corresponding effective temperature is (roughly) 80,000 K, 27,000 K, and 18,000 K, 
respectively; the models thus cover a representative range of effective temperatures of the AGN 
disk annuli. 

For simplicity, we consider disks composed of hydrogen and helium only. This allows us to 
study various effects without spending too much computer time, while taking into account all 
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the essential physics. The effects of heavy elements on disk models will be considered in the 
subsequent paper. Here, we intend to investigate two questions, namely, i) the effect of degree of 
sophistication in the disk modeling (LTE versus NLTE; effects of lines); and ii) the sensitivity of 
computed model structure on the individual viscosity parameters. 

Hydrogen is represented essentially exactly: The first 8 levels are treated separately, while 
the upper levels are merged into the averaged non-LTE level accounting for level dissolution as 
described by Hubeny, Hummer, & Lanz (1994). Neutral helium is represented by a 14-level model 
atom, which incorporates all singlet and triplet levels up to n = 8. The 5 lowest levels are included 
individually; singlet and triplet levels are grouped separately from n = 3 to n = 5, and we have 
formed three superlevels for n = 6,7, and 8. The first 14 levels of He"*" are explicitly treated. We 
assume a solar helium abundance, A^(He)/A^(H) = 0.1. 

In the NLTE/L models, all the lines are treated explicitly, assuming Doppler profiles with 
turbulent velocity ranging from the thermal velocity up to the vertically-averaged sound speed, 

Vturh = \/^, (48) 

where (-P/p) is given by Eq. (|l6|). For the three annuli at r = 2, 11, and 20, of the disk specified 
above, the sound speed has values of 3880, 770, and 375 km s~^, respectively. Although considering 
different values of Vturb has a significant effect on computed rest-frame line profiles, we found a 
negligible effect on the resulting disk structure. 

3.1. Non-LTE Effects 

Figure 1 shows the temperature as a function of depth for the three annuli. The depth is 
expressed as column mass in g cm^^. We see several interesting features. Firstly, for the hot 
model at r = 2 (Tgfj = 80,000 K), the NLTE temperature structure differs appreciably from the 
LTE structure, even at the midplane. The cooler models exhibit different temperature at the 
upper layers, i.e. for logm ~ 1.5, while the temperature structure is unchanged by NLTE effects 
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in the deep layers. This is explained by the fact that the effective optical thickness 

TeS ~ V'^tot Tabs , 

(where rtot is the total optical depth corresponding to %, i.e., including electron scattering, while 
Tabs is the optical depth corresponding to the true absorption, k), becomes comparable to or 
smaller than 1 for hot models. In other words, a photon created at the midplane has a non-zero 
probability that it will undergo a series of consecutive scatterings without being destroyed by a 
thermal process until it escapes from the disk surface. Consequently, even the deep layers of the 
disk now effectively "feel" the presence of the boundary, which gives rise to NLTE effects even 
close to the midplane. 

Secondly, as expected, differences between NLTE/C and NLTE/L models are negligible in 
the deep layers, while they are important in the upper layers - the lines heat up the upper disk 
atmosphere (for log m ~ 0.5) by some 10,000 K for the hot model; by 4,000 K for the intermediate 
(r = 11) model, and by only 1000 K for the "cool" model. The effect is exactly analogous to the 
temperature rise predicted in hot stellar atmospheres (Mihalas 1978), which is explained as an 
indirect effect of Lyman and Balmer lines on the heating in the Lyman and Balmer continuum. 
The explanation goes as follows: roughly speaking, the temperature structure is determined by 
the balance between the radiative heating, which is mostly provided by the Lyman (for hotter 
models) and the Balmer (for cooler models) continuum, and the radiation cooling in the free-free 
transitions in the optical and infrared region. Radiative transfer in the Lyman (Balmer) lines 
explicitly gives rise to an overpopulation of the hydrogen n = 1 (n = 2) states, which leads to an 
increase of the efficiency of the Lyman (Balmer) continua which in turn leads to an additional 
heating at the surface. There is a competition between this heating and the traditional surface 
cooling caused by the lines, but in the present case the indirect effect dominates. 

In Fig. 2, we plot the density structure for the same models as displayed in Fig. 1. For the 
hot model, density is almost unaffected by NLTE effects. This is easy to understand, since in this 
case the radiation pressure is completely dominant, and the disk structure is given by the following 
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simple analytic formulae. The hydrostatic equilibrium equation ( pSD reduces to 



dm 



Qz, (49) 



where we have denoted Q = {GM/R ) {C /B). Using Eq. (35) and the first moment of the transfer 



equation, (29), we may express the gradient of the radiation pressure as 



du = — / — Hydv = — xh tl , (50) 



dm c Jo dm c Jo p c 

where xh is the flux-mean opacity, defined by 



oo 



XH= / {Xu/p)H,dv/H , (51) 





and H is the total (frequency- integrated) Eddington flux, H = H^du = AttFj.^. The 
flux-mean opacity is particularly simple in the case when the electron scattering is the 
dominant source of opacity. In this case xh = neae/p- For a medium where hydrogen 
is fully ionized and helium is partially doubly ionized, the flux-mean opacity is given by 
XH = (o"e/m//)(l + aY)/{l + 4y) « 0.4 (1 + aY)/{l + iY), where Y is the helium abundance, 
Y = A^(He)/A^(H), and a = 1 + A^(He''~''')/A^(He); mn is the mass of hydrogen atom. 

Using Eq. (^), we obtain 

^T^^[l-9{m)] = Qz. (52) 
c 

We know that at the disk surface, m = 0, 6{m) = 0, and therefore the total disk height, h, is given 
by 

The hydrostatic equilibrium equation thus reduces to 

^ = l-e{m) (54) 

This equation enables us to estimate the density, since from Eq. ( |4^ ) we know that 1/p = —dz/dm, 
so that from Eq. (jsj) we obtain 1/p = h{d9/dm). However, from the definition of 6, Eq. (pO|) we 
have d9/dm = w{m) / iwmo), so that we finally obtain 

p{m) = (55) 
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This equation shows that the density structure for the case of dominant radiation pressure is solely 
determined by the dependence of viscosity on depth. For a depth-independent viscosity, density is 
constant, and is given by p = mo/h, i.e the total column mass divided by the disk height. 

Before proceeding further, we derive the term {P/j)) for the case of radiation-pressure 
dominated disks. Assuming a depth-independent viscosity, the density is also constant, 
p = 'p. The hydrostatic equilibrium equation thus reads dP/dz = —Q'pz, which has solution 
P{z) = Q'p{h'^ — z^)/2, and consequently P = Q'ph'^/3. Substituting for h from Eq. (|53|), we 
finally obtain 

P ^ 3GMM f XH y , . 

p \87rmHcJ BC ^ ' 

Substituting the pure-hydrogen form of XH, XH = Ce/""!-// into Eq. (p6|), we obtain Eq. ( p^ ) used 

in the definition of uq. 



It should be stressed that Eqs. through (55) apply if the gradient of the radiation pressure 



dominates over the gradient of the gas pressure, dP^a^^^/dm » dPgs^/dm; it is not sufficient that 
^rad ^ Pgas- It is clear that even in the case where the radiation pressure dominates everywhere 
in the atmosphere, the radiation pressure gradient becomes very small in the upper layers where 
the optical depth is smaller than unity, because the radiation field is already formed and does 
not change when going outward. The gas pressure gradient then takes over in these layers, and 
consequently density starts to decrease exponentially with height. For a comprehensive discussion, 
see Hubeny (1990a). 



In any case, Eq. (55) shows that in the inner layers of the radiation-pressure-dominated 
disks, density structure does depend only on viscosity, and is therefore insensitive to NLTE 
effects. This is indeed demonstrated in Fig. 2. The above considerations also explain that for 
the radiation-pressure-dominated disks the disk height depends only on the radial distance R, 
through the relativistic correction D/C, and slightly through XH (because of variations in the 
ionization of helium). The z-m relation is given through the depth dependence of viscosity, and is 
thus insensitive to local temperature, as well as to NLTE effects. This is demonstrated in Fig. 3, 
where we plot the vertical distance z as a function of m for the same models as displayed in Fig. 1. 
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Figure 4a displays the departure coefficient (6-factor) for the ground state of hydrogen for the 
three annuli, and for NLTE/C and NLTE/L models. Departure coefficient is defined as 6j = rij/n*. 
Near the midplane of the disk, the 6-factor is close to unity because the optical depth is large. It 
starts to deviate from unity as soon as the Lyman continuum becomes effectively optically thin. 
The behavior of the hot model differs from the two cooler ones. 

In the hot model, electron scattering dominates the opacity even in the Lyman continuum; 
the thermal coupling parameter 



is very small in the inner layers (e.g., « 2 x 10~^ at the midplane). Consequently, starts to 
deviate from the thermal source function, which is roughly equal to the Planck function, already 
in deep layers - see the upper left panel of Fig. 7. One can make these considerations more 
quantitative by invoking a simple model of radiative transfer in the presence of scattering (see, 
e.g., Mihalas 1978). In a simple case of depth- independent B^, and e^,, the mean intensity is given 



It is clear that By for ~ i/Se^ (called the thermalization depth) ; while for < \f^e^ the 

mean intensity Jy drops below By\ the drop is larger for smaller ey. We note that at the surface, 

Ji/(0) ~ By. 

The photoionization rate is proportional to an integral of J,y, while the recombination rate 
is proportional to an integral of By. As a result, the recombinations dominate over ionizations, 
which leads to an overpopulation of the hydrogen ground state. For cooler models, the formation 
of Lyman continuum is different. The electron scattering is not overwhelmingly dominant, so 
the thermal coupling parameter e is now smaller (e is between 10~^ and 10~^ in the midplane 
layers), and the mean intensity decouples from By farther away from the midplane. The 
Planck function at the frequencies of the Lyman continuum decreases very fast with decreasing 
temperature, because the Lyman continuum frequencies are in the Wien tail of the Planck function, 
By{T) oc exp{—hi>/kT). This decrease in the local value of By is now faster than the decrease 



by 



Jv{Tu) ^ By{Ty) 1 + ^ - ex.^{-\/^Ty) / {l + ^/€^) . 



(58) 
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of Jjj which, when going outward, is more and more decoupled from the local temperature. We 
thus obtain Jy > By throughout most of the atmosphere; consequently the ionizations dominate 
over recombinations in the Lyman continuum, and hence the hydrogen ground state becomes 
underpopulated . 

Once the Lyman continuum becomes transparent (which happens around logm ^ for all 
models), hi does not change much for NLTE/C models. For NLTE/L models, however, 61 starts 
to increase outward above logm « 0, which is caused by radiative transfer in La and other 
Lyman lines. This is a standard NLTE effect - strong resonance lines cause the lower level to be 
overpopulated, while the upper level becomes underpopulated with respect to LTE (e.g., Mihalas 
1978). This is illustrated in Fig. 4b, which displays the 6-factor for the n = 2 state of hydrogen. 
Detailed explanation of the behavior of 62 is quite analogous to that of 61 discussed above. 

The most interesting quantity predicted by the models is the emergent radiation. Figure 5 
compares the emergent flux in the Lyman limit region for the three models, LTE, NLTE/C, and 
NLTE/L. Only hydrogen and helium lines are taken into account. In the NLTE/C model, the 
line source functions are computed using NLTE level populations, which in turn were calculated 
by TLUSDISK by considering the lines to be in detailed radiative balance. Consequently, the 
predicted line profiles are inconsistent, and we show them only for demonstration purposes. 
The spectra were computed using program SYNSPEC (Hubeny, Lanz, k. Jeffery 1994), which 
calculates synthetic spectra for model atmospheres or disks previously computed by TLUSTY or 
TLUSDISK. We assume the turbulent velocity Uturb = 3880, 770, and 375 km s~^, for the annuli 
at r = 2, 11, and 20, respectively. For a clearer display, the final spectra are convolved with a 
gaussian with FWHM = 10 A. We stress that the spectra are computed in the rest frame of the 
annulus; to obtain spectrum received by a distant observer one would have to take into account 
Dopplcr velocities, and general rclativistic effects (frequency shift, light bending, etc.). As pointed 
out in the Introduction, we will study here only the rest-frame radiation. 

For the hot model, the NLTE/C and NLTE/L models give essentially the same flux across the 
Lyman discontinuity which in this case virtually disappears. In contrast, the LTE model predicts 
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the Lyman jump in emission. In the intermediate model, LTE predicts almost non-existent 
Lyman jump while NLTE/L predicts a weak emission jump. The NLTE/C model predicts a 
larger emission jump, and spuriously strong emission in the Lyman lines. In the cool model, LTE 
predicts the Lyman jump in very strong absorption, while the NLTE/L models again predicts 
essentially smooth spectrum in this region. 

The explanation of this behavior follows from the above discussion of the hydrogen departure 
coefficients bi and 62 (Fig. 4), and from studying the behavior of By, Sy, J^, and as functions 
of depth. This is displayed in Fig. 7 for the hot model at r = 2, and in Fig. 8 for the cool model 
at r = 20. Let us ffi'st consider LTE models. In the hot model, electron scattering completely 
dominates, so the total opacities at both sides of the Lyman discontinuity are roughly equal; 
hence both sides of the discontinuity are formed at the same depth. The only difference is the 
value of the thermal coupling parameter e^. Since the thermal opacity on the blue side (f;,) of the 
discontinuity is larger that that on the red side (i^r), we have e^^ > e^^. Consequently, the mean 
intensity at the red side is more uncoupled from By, and consequently Jy^^ > Jy^ - see Eq. (^Sj). In 
fact, this is an interesting manifestation of the classical Schuster mechanism. By the same token, 
one can understand why the Lyman lines also appear in emission in the hot LTE model. 

In the cool model, in contrast, we have a classical LTE jump: the high-opacity (blue) side of 
the jump is formed much higher than the low-opacity (red) side, because the electron scattering 
opacity no longer dominates over the thermal opacity. Since the temperature decreases outward, 
the blue side of the jump is formed at lower temperature, and the flux is consequently lower. In 
the intermediate model, both effects compete, but the first one wins, so we obtain a weak emission 
jump. 

NLTE effects modify this picture significantly. In the hot model, the red side of the Lyman 
jump remains almost unchanged by NLTE effects because the departure coefficient 62 is close to 
unity around the thermalization depth of the Balmer continuum. In the Lyman continuum, the 
thermal coupling parameter ey^^ > Ey^ as discussed above. However, in NLTE one has to modify 
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Eq. (ID to read 

where is the thermal source function. In the case of Lyman continuum, it is roughly given 
by ~ By/bi. Since 61 > 1 in the continuum-forming layers (see Fig. 4a), we have S'^ < By. 
We see that the effect of an increased thermal coupling parameter at the Lyman continuum 
frequencies (which tends to increase the emergent intensity) is offset by a decrease of the thermal 
source function. The latter results from a predominance of recombinations over ionizations, as 
explained above. 

Finally, we explain why we obtain a somewhat higher flux in the Lyman continuum 
for the NLTE/C models than for NLTE/L models. This is given by the fact that 
6i(NLTE/C) < 61 (NLTE/L) all the way from the surface to deep layers (see Fig. 4a), 
which in turn is the effect of La and other Lyman lines. Because of low density, the photon 
destruction parameter e is rather small and therefore the thermalization depth is rather large - 
much larger than the depth of formation (see, e.g., Mihalas 1978). Consequently, the thermal 
source function in the Lyman continuum, ~ B^/bi, is smaller in the NLTE/L model, and 
consequently the emergent flux is lower. 

Figure 6 displays the EUV continuum. As we have discussed in Paper I, NLTE models 
produce the He II Lyman jump in emission for the hot annuli. For cooler annuli, both LTE and 
NLTE models produce a strong absorption in the He II Lyman jump. Both NLTE/C and NLTE/L 
models produce very similar results. 

Having computed models of vertical structure, we can verify that neglecting Compton 
scattering was indeed a legitimate approximation. To this end, we compute the appropriate 
Compton y parameter which specifies whether a photon will be significantly influenced by 
comptonization in traversing the medium. For non-relativistic electrons, the Compton y parameter 
is given by (Rybicki & Lightman 1979) 

y = max (res, r|,) (60) 



1 + - exp(- 



(59) 



- 26 - 



where t^s is the total electron-scattering optical depth of the medium. In the case in which thermal 
absorption is not negligible, the t^s depends on frequency, and is defined by (Rybicki & Lightman 
1979, their eq. 7.42) 

resile) - (tT^) = (1 - ' (61) 

where the second equality follows from Eq. (|57|). Strictly speaking, the above formulae are derived 
for a homogeneous medium. We may nevertheless use them for some characteristic (averaged) 
values of structural parameters. The Compton scattering would be most important in the hottest 
model. Taking T = Tgfj and the characteristic e for the Lyman continuum frequencies equal to 
10^^ (see Fig. 7), we obtain for y ~ 5 x 10^'^. Even the upper limit for y, taking the midplane 
values of r ~ 3 X 10^ and e ~ 2 x 10^^ (which would give a largely overestimated value of y 
because the local temperature is lower and the parameter e is higher) yields y ~ 10^^, which is still 
well below unity. Therefore Compton scattering is never very important for the models considered 
here. 

In conclusion, there are several NLTE effects; for hot and cool models they decrease the 
magnitude of the Lyman discontinuity as compared to LTE models, or to simplified NLTE/C 
models, while for the intermediate temperatures they somewhat increse the Lyman jump. We also 
find that for hot annuli, the NLTE/C models are sufficient for predicting the emergent continuum 
radiation, while for cooler annuli the NLTE/C models produce a somewhat higher flux in the 
Lyman continuum, and spuriously strong emissions in the Lyman lines. 

We stress again that the present paper is devoted to studying a few representative annuli of an 
AGN accretion disk model. We do not aim here to answer the important question of whether the 
present models will alleviate the long-standing Lyman jump problem, namely that the theoretical 
models predict a significant jump, while the observed jump is virtually non-existent. The problem 
arose from early calculations by Kolykhalov & Sunyaev (1984); a comprehensive review of the 
current status was presented by Blaes (1998). To address this problem, we have to compute 
theoretical models and emergent spectra for all annuli of a disk and integrate them using the 
appropriate general relativistic transfer function (Cunningham 1975; Agol 1997). Also, we should 
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take into account effects of metal lines, which was recently found to be quite important (Hubeny 
& Hubeny 1998; Agol, Blaes, & Hubeny 1998). We shall defer this study to a future paper of this 
series. 

3.2. Effects of Viscosity 

The most important parameter influencing the model structure is oq. We first study the 
effects of changing oq, while keeping the remaining viscosity parameters fixed. We take Co = 
(i.e., a constant viscosity in the inner region); = 10^^ (i-e., the viscosity starts to decrease 
with m only for m/rriQ = 10~^; in other words, 99% of the mass of the disk column has a constant 
viscosity); and (i = 2/3. Since the total column mass, mo, is proportional to 1/ao - see Eq. ([T^), 
we may think of effects of changing the value of the parameter ao being in fact effects of changing 
the disk column mass. We stress again that our study concerns the behavior of the individual disk 
annuli and its sensitivity to the adopted value of oq. In reality, it is not clear whether the same 
value of ao should be considered for all annuli of a given disk or not. This will only be solved by 
future detailed MHD simulations similar to those of Stone et al. (1996). 

Figure 9 shows the results for the hot annulus; the upper panel displays the temperature, 
the middle panel the density, and the lower panel the emergent flux, for the NLTE/C models for 
ao = 0.3, 0.1, 0.03, 0.01, and 0.005. The behavior of density is easy to understand. The total 
column mass, uiq, is directly proportional to 1/ao ~ see Eq. (|l8|). The total disk height does not 
depend on viscosity for radiation-pressure-dominated disks (see Eq. ^3|); therefore the density in 
the inner regions is proportional to niQ - see Eq. ( |55|) - and consequently p{mK,mQ) oc 1/ao- This 
is demonstrated on the middle panel of Fig. 9. Consequently, the thermal coupling parameter e 
(which is proportional to density), is larger for smaller uq. 

The behavior of temperature is more complicated, since it follows from several competing 
mechanisms. First, we write a simple analytic formula derived by Hubeny (1990a) which gives a 
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reasonably accurate estimate of the local temperature, viz., 

1 inf m) 

(62) 



1 1 wim) 



2rtot/ VS 3 mo Kb (t) w 
where r is the flux-mean optical depth, dr = xudz = (xH / p)dm; rtot is the flux-mean optical 
thickness at the midplane; and kb{t) is the Planck-mean opacity, defined by 

KB = {Ky/p)Bydu/B , (63) 



where B is the frequency-integrated Planck function, B{T) = Bu{T) du = {(t/tt)T'^. This 
formula is easy to understand. The term r -|- \/V?> in the square bracket is the same as for the 
classical LTE-grey semi- infinite model stellar atmospheres in radiative equilibrium (i.e., no energy 
generated in the atmosphere). The "correction" (1 — r/2rtot) reflects the fact that the disk has a 
finite total optical thickness; while the third term describes the energy generation in the disk due 
to viscous dissipation. Notice also that the Planck-mean opacity contains the thermal absorption 
coefficient, Kj/, while the flux-mean opacity contains the total absorption coefficient, Xu- The fact 
that the dissipation term of Eq. ( |62|) contains the Planck mean opacity is easy to understand, 
since it is the thermal absorption, not scattering, that contributes to the global energy balance. In 
the case of dominant electron scattering, the relation between the flux-mean optical depth, r, and 
the column mass, m, is particularly simple, r ~ 0.34m (if helium is doubly ionized), or r ~ 0.31m 
(if helium is singly ionized), as discussed in Sect. 

As follows from Eq. (B^), the temperature at the midplane is roughly given by 



Ttot ^ 1 



(64) 



2 3 mo KB (rtot) 

where we assume that the disk is optically thick, rtot ^ 1, and the local viscosity at the midplane 
roughly equal to the averaged viscosity. Quantity mo kb (rtot) is very roughly equal to the total 
Planck-mean optical depth of the disk. We write moKB(rtot) = ertot. Since e is proportional to 
k,b/xHi it may be called the "averaged thermal coupling parameter". The midplane temperature 
then becomes 

TV«) = |^»(^ + 5^). (05) 
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If the electron scattering is negligible, e ~ 1, and the second term in Eq. I pq ) is negligible compared 
to the first one in the case of optically thick disks. In fact, the disk behaves like a normal stellar 
atmosphere. However, as soon as the averaged thermal coupling parameter becomes small, the 
second term in Eq. (|6^) starts to contribute. This is indeed illustrated in the upper panel of 
Fig. 9. For the lowest ao, cto = 0.005, the total optical depth is very large, and consequently the 
central temperature is high. Increasing oq, the total optical depth decreases, and so does the 
central temperature. However, with increasing ao the density decreases, and therefore the thermal 
coupling parameter e also decreases. Hence the second (viscous-heating) term of Eq. ( |65| ) starts to 
make appreciable contribution to the central temperature, which first levels off and then starts to 
increase even if the total optical depth decreases (see the last curve for ao = 0.3). Since the second 
term of Eq. (65) makes roughly the same contribution at all depths, we see a similar increase of 



local temperature at all depths. 

At the top of the disk, density decreases exponentially, and consequently e decreases faster 
than w{m)/W, which is assumed to be a simple power law (recall that w{m)/w oc (m/mo)^/^ for 
the models displayed in Fig. 9). Since for the model with the largest ao, ao = 0.3, the averaged 
thermal coupling parameters is so small that the second term in fact determines the temperature 
structure. Consequently, the surface temperature begins to increase to large values (this effect 
is better seen in Fig. 12 where we consider a lower (^i and therefore a higher viscosity in the 
upper layers). However, as discussed by Hubeny (1990a), this increase is partly spurious, because 
in computing the Planck-mean opacity we have included only lines (and continua) of hydrogen 
and helium. In reality, however, there is a host of other lines of light elements which operate in 
the upper layers (e.g., the resonance lines of C IV, N V, O VI, and higher ions of S, Ne, Fe, 
etc., for even higher temperatures). Including these lines will increase the Planck-mean opacity 
significantly, which physically means that we are including efficient mechanisms which are able to 
radiate the dissipated energy away. 

The lower panel of Fig. 9 shows the continuum flux for all models. The Lyman jump varies 
from being almost non-existent for the low-ao models, to a weak absorption jump for the high-ao 
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models. As explained above, the appearance of the Lyman jump is a result of a competition of 
two mechanisms: the overpopulation of the ground state of hydrogen, which decreases the source 
function and therefore the emergent flux in the Lyman continuum; and a higher thermal coupling 
parameter which causes the mean intensity in the Lyman continuum frequencies to be more 
coupled to the thermal source function, and therefore causes the emergent flux to be higher. In 
the present case, going to higher oq decreases density in the inner layers, therefore causing the 
hydrogen ground state to be more and more overpopulated, and since this is a dominant effect, 
the flux in the Lyman continuum decreases with increasing oq, so the Lyman jump is driven to a 
somewhat stronger absorption. 

In contrast, the magnitude of the He II Lyman jump decreases with increasing oq, from being 
a very conspicuous emission jump at ao ~ 0.03 to a modest emission jump at oq = 0.3. This is 
explained by the fact that the strong jump at ao ~ 0.03 is essentially caused by the Schuster 
mechanism discussed above. When going to higher ao the density decreases, and the electron 
scattering is more and more important; the jump thus becomes weaker. However, it is important 
to realize that the number of He II ionizing photons remains roughly the same for all models, since 
the flux at high frequencies {u ~ 1.7 x 10^^ Hz) is higher for higher ao; this is a consequence of 
higher temperature in the continuum-forming region - see the upper panel of Fig. 9. 

Figure 10 presents an analogous comparison of models with various values of ao for a 
"cool" annulus, Tgff = 18, 000 K. The behavior of models is similar to that of the hot annulus, 
although the effects of viscosity are generally weaker. The total column mass and the inner 
density are proportional to 1/ao, as in the previous case. The central temperature now decreases 
monotonically with increasing ao because the electron scattering is not dominant and therefore 
the second term in Eq. (|65| ) is not very important. The Lyman jump varies from a weak emission 
for high-ao models to a weak absorption for low-ao models. This follows from the fact that for 
increasing ao the density decreases, so the magnitude of NLTE effects increases. In the present 
case, the ground state of hydrogen becomes more underpopulated, and consequently the flux in 
the Lyman continuum increases. 
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Next, we examine the effects of changing the other viscosity parameters. In Fig. 11, we display 
the models of the hot annulus computed for fixed oq (ao = 0.1) and the power law exponents 
(Co = and ("i = 2/3), and for several values of the division mass, ma. The division mass is varied 
from rn-d = 1 (i.e., no inner region of constant viscosity), to = 0.3, 0.1, 0.03, and 0.01. As 
expected, the only interesting effect is the behavior of density in the inner layers, which is governed 
by Eq. (^5|). The central temperature is the largest for = 1 because the central density, and 
thus the thermal coupling parameter, are lowest. The emergent flux is only weakly influenced by 
changing mj; a modest change of the flux is seen in the hydrogen and He II Lyman continuum 
which are the most sensitive to the thermal coupling parameter. 

Finally, Fig. 12 displays the effect of changing Ci, Ci = 2/3, 1/2, and 1/3, for oq = 0.1, 
nid = 1, for the hot annulus (r = 2). For Ci ~ 1/2, the only appreciable effects are seen in the 
behavior of density in inner layers. The overall behavior of the models is easily explained by the 
same considerations as above. For the model with the lowest Ci; Ci = 1/3, the local viscosity 
in the upper layers is so high that we see the effects of temperature runaway at m ~ 5 x 10"'^. 
Decreasing the value of ("i still would increase this instability. However, as discussed above, we 
cannot study this effect with simple H-He models, and we leave this problem to a future paper. 



As discussed in Sect. 2.2, we consider the depth-dependent kinematic viscosity as a step-wise 
power law function given by Eqs. ( p2|) - (|25|) . In contrast, the local a-prescriptions, e.g., the 
variant suggested by Dorrer et al. (1996), considers the kinematic viscosity in the form given by 



This approach essentially considers the viscosity being proportional to the gas pressure at low 
optical depths, while being proportional to the radiation pressure at large depths; r is the 
Rosseland mean optical depth. 



3.3. Comparison with the Local a— viscosity Approach 



Eqs. (p!^), with the turbulent velocity given by 




(66) 
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We can compute an "effective a-parameter" in such a way that the two prescriptions give the 
same value of local viscosity, 

aofr(m) = — — . (67) 

n-i'turbim) 

Since our parameterization does not take viscosity to be proportional to the local turbulent 
velocity, the resulting agff will be depth-dependent. 

In Fig. 13, we plot Oeff for the hot annulus (r = 2), and for various viscosity parameters. We 
see that the values of Oeff for our standard model, ao = 0.1, Co = 0, Ci = 2/3, and ra^jra^ = 0.01 
are located in a reasonable range of 0.02 ~ Oeff ~ 0.5. The behavior of Ocs as a function of depth 
is easily understood. In the inner layers, Uturb = {P/pY^"^ ^ (-frad/p)^^^- Since the optical depth 
is large, the radiation pressure can be approximated by the thermodynamic equilibrium form, 
Prad oc . In the model with rrid/mQ = 0.01 density is roughly constant for nid ~ m < ttiq. 
Consequently, Uturb oc T^. To first order, oc t oc m (see Eq. 62), so that finally Vturb m^/^. 



and Oeff oc m^^/^ for m > m^. This is indeed seen in the upper and lower panels of Fig. 13; Oeff 
decreases with m somewhat slower than m^^/^ because the temperature increases slower than 
oc m - in fact, oc (m — m^/2mo). 

For models with a single power law viscosity (i.e., = mo), the behavior of a^s in the 
inner layers can also be easily understood. Here w{m) oc mf' (we write C, instead of C,i to simplify 
the notation), so the density varies as p{m) oc l/w{m) oc m~^. The radiation pressure scales, as 
discussed above, as Prad oc m. Consequently, vturb m*^^+^)/^, and thus Oeff oc ■m^^~^^^'^ . This is 
demonstrated in the middle panel of Fig. 13; for instance, we see that for C = 1, a^s is indeed 
almost constant in the inner layers. In the outer layers, the density decreases outward faster than 
w{m), so OeS decreases faster than m^'^^^^^'^. 

In the gas-pressure-dominated regions, fturb = (-fgas/p)^^^i so we have fturb oc T, and thus 
Oeff oc w{m) /T{m). In the outer layers, the temperature is roughly constant with m, so that 
Oeff oc w{m) oc m''i. Again, this is clearly seen in Fig. 13. 

The important point to realize is that although the local value of varies significantly, its 
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infiuence on the disk structure is rather small. Consider for instance the middle panel of Fig. 13, 
and compare it to Fig. 12, which displays the same models. Although cteff differs by several orders 
of magnitude, the disk structure is hardly affected. From the point of view of constructing detailed 
models of vertical structure of AGN disks this is a good news, because the most uncertain part of 
physics - the viscous dissipation - has relatively small effect on the computed structure. However, 
we should bear in mind that this study was limited to considering a "well-behaved" viscosity which 
smoothly decreases towards the disk surface. When one assumes, for instance, that the viscous 
dissipation is concentrated mostly in the outer layers (e.g., Sinccll &: Krolik 1997), the overall 
vertical structure may be significantly different. We plan to study such non-standard models in 
future papers of this series. 

4. CONCLUSIONS 

We have calculated several representative models of vertical structure of an accretion disk 
around a super massive Kerr black hole. The interaction of radiation and matter is treated 
self-consistently, taking into account departures from LTE for calculating both the disk structure 
and the radiation field. The viscosity is parameterized through the parameter oq that describes 
the vertically averaged viscous stress, and two power-law exponents and ("i, and the division 
point md between these two forms. The disk structure and emergent radiation is sensitive mainly 
to the values of ckq, while the other parameters influence the disk structure to a much lesser extent. 
However, although the detailed shape of the predicted spectrum is sensitive to adopted ckq, the 
overall appearance of the spectrum is quite similar in all cases. 

We have shown that effects of departures from LTE are very important for determining the 
disk structure and emergent radiation, particularly for hot and electron-scattering dominated 
disks. We have shown that at least for the disk parameters studied in this paper, NLTE effects 
typically tend to diminish the value of the Lyman jump; in hot models they suppress the Schuster 
mechanism by which the LTE models produce a strong emission jump, and in cooler models they 
increase the flux in the Lyman continuum due to an underpopulation of the hydrogen ground 
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state. Also, we have shown that relaxing the approximation of detailed radiative balance in the 
hydrogen and helium lines (i.e., computing the so-called NLTE/L models) changes the predicted 
line profiles significantly, but otherwise does not yield significant changes in computed vertical 
structure or emergent continuum flux. This result shows that for estimating the continuum 
radiation of AGN disks composed of hydrogen and helium, the NLTE/C models provide a 
satisfactory approximation. 

So far, we have limited our analysis to a simple H-He chemical composition. A preliminary 
study (Hubeny & Hubeny 1998) indicates that the effects of numerous metal lines on the prediced 
spectral energy distribution of AGN disks may be quite significant. However, that study used 
a simplified approach in which the vertical structure was fixed by a H-He model, while the line 
opacity was taken into account only in the spectrum synthesis, assuming LTE source function 
in metal lines. Such an approach is inconsistent, since, first, the metal lines may change the 
disk vertical structure (the so-called metal line blanketing effects, long known form the theory 
of classical stellar atmospheres) and, second, the source function in metal lines may depart 
significantly from the LTE value. We will therefore need to construct self-consistent, fully 
metal-linc-blankctcd models of vertical structure of AGN disks, taking into account effects of 
literally millions of spectral lines in NLTE. A work of this project is under way, and will be 
reported in a future paper of this series. 

The results presented here do not indicate any fatal flaw of the AGN accretion disk paradigm. 
In contrast, they show that one of the previous critical arguments against the accretion disk 
paradigm, the magnitude of the Lyman jump, essentially disappears when increasing a degree of 
realism of the modeling procedure by relaxing previous simplifying approximations, in particular 
the local thermodynamic equilibrium and a simplified vertical disk structure. However, this 
study has concentrated on only one aspect of the problem, the spectral energy distribution in the 
optical, UV, and EUV region. Many questions, such as the overall spectral energy distribution of 
the whole disk, the effects of external irradiation, the continuum polarization, etc., remain to be 
explored in detail. This is exactly what we intend to do in future papers of this series. 
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Fig. 1. — Temperature as a function of depth for three annuh at r = 2 - upper panel; r = 11 
- middle panel; and r = 20 - lower panel; for a disk model with the mass of the black hole, 
M = 2 X 10^ Mq, the mass accretion rate M = IMq/jt, and the maximum stable rotation, 
a/M = 0.998. The viscosity parameters are taken ckq = 0.1, (q = 1, d = 2/3, and = 0.01. 
For all annuh, the thick line is the NLTE/L model, the dashed line the NLTE/C (i.e., NLTE with 
continua only) model, and the dotted line is the LTE model. 
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Fig. 2. — Mass density (in g cm~^) for the model annulus at r = 2 - the upper curves, and for 
r = 20 - the lower curves; for the same models as displayed in Fig. 1. We did not show the r = 11 
model because its behavior is quite analogous to that of the r = 20 model. 
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Fig. 3. — Vertical distance from the disk plane as a function of column mass for the model disk 

annulus at r = 2 (lower curves); r = 11 (middle curves); and r = 20 (upper curves), for the same 
models as displayed in Fig. 1. Notice that NLTE effects upon the z vs. m relation are quite small. 
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Fig. 4. — NLTE departure coefficients for tfie hydrogen ground level (a), and tfie n = 2 level 
(6); for the NLTE/C models (dotted lines), and NLTE/L models (full lines). The lines without 
additional symbols correspond to the model annulus at r = 2; the lines with additional "+" signs 
to the r = 11 models, and with diamonds to r = 20 models. 
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Fig. 5. — A comparison of emergent flux in the region around the Lyman discontinuity for the 
NLTE/L model (full line); NLTE/C model (dashed line), and LTE model (dotted line) of the model 
disk annulus at r = 2 (upper panel), r = 11 (middle panel), and r = 20 (lower panel). 
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Fig. 6. — The same as in Fig. 5 for the EUV spectrum at the region of He II Lyman discontinuity 
(A = 227 A), and the He I ground-state discontinuity (A = 504 A). 
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Fig. 7. — Plot of the mean intensity of radiation (full line, labeled J), the Planck function (dashed 
line, labeled B), the thermal source function (dot-dashed line, labeled S), and the thermal coupling 
parameter e (dotted line), as a function of the monochromatic optical depth, for the model annulus 
at r = 2. The two upper panels are NLTE models; the two lower panels are LTE models. The left 
panels display the values for the blue side of the Lyman discontinuity (A = 900 A), and the right 
panels display the values for the red side of the Lyman discontinuity (A = 950 A). 
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Fig. 8. — The same as in Fig. 7, but for the model annulus at r = 20. 
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Fig. 9. — The effect of oq on the vertical structure and emergent flux for the model annulus at 
r = 2. Other basic parameters are the same as for the models displayed in Fig. 1, namely Co = 1; 
Ci = 1, and rrij = 0.01. Upper panel: temperature; middle panel: mass density; and lower panel: 
emergent flux. The thick line corresponds to the highest ao, ao = 0-3, the line with additional "+" 
signs to the lowest one, ao = 0.005, and the curves in between to ao = 0.1, 0.03, and 0.01. The 
curves for the two extreme values are labeled by the values of ao- 
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Fig. 10. — Analogous to Fig. 7, for the model annulus at r = 20; the adopted ao values are 0.45 
(thick line), 0.3, 0.1, 0.03, and 0.01 (the line with additional "+" signs). 
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Fig. 11. — The effect of (the division point between two different power laws of viscosity 
parameterization) on the vertical structure and emergent flux for the model annulus at r = 2. 
Other basic parameters are the same as for the models displayed in Fig. 1, namely ("o = 1) Ci = 2/3, 
and ao = 0.1. The panels are arranged as in Fig. 7. The adopted values of ma are: 1 (thick line); 
0.3; 0.1; 0.03; and 0.01 (the line with additional "+" signs). The curves for the two extreme values 
are labeled by the values of mj. 
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Fig. 12. — The effect of the power law exponent Ci on the vertical structure and emergent flux 
for the model annulus at r = 2. Other basic parameters are ckq = 0.1, and = 1 (i.e., no inner 
region of constant viscosity). The panels are arranged as in Fig. 7. The adopted values of Ci are: 1 
(thick Hne); 2/3; 1/2; and 1/3 (the line with additional "+" signs). The curves for the two extreme 
values are labeled by the values of . 
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Fig. 13. — Plot of effective a-viscosity parameter that would yield the same value of kinematic 
viscosity in our models and in the a-parameterization of viscosity suggested by Dorrer et al. (1996). 
Upper panel - Oeff for the models displayed in Fig. 7, i.e., the models of the r = 2 annulus with 
Co = 0, Ci = 2/3, and nid = 0.01, and with various values of ao- The adopted values of aQ are (from 
top to bottom): 0.3 (thick line), 0.1, 0.03 and 0.01 (the line with additional "+" signs). Middle 
panel - a^s for the models displayed in Fig. 10, i.e., the models of the r = 2 annulus with ckq = 0.1, 
rrid = 1, and with various values of Ci- The adopted values of (i are (from bottom to top): 1 (thick 
line); 2/3; 1/2; and 1/3 (the Hne with additional "+" signs). Lower panel - Ocff for the models 
displayed in Fig. 9, i.e., the models of the r = 2 annulus with = 0.1, Co = 0, Ci = 2/3, and with 
various values of the division depth md- The adopted values of rrid are (from bottom to top): 1 
(thick line); 0.3; 0.1; 0.03; and 0.01 (the line with additional "+" signs). 






